home *** CD-ROM | disk | FTP | other *** search
/ io Programmo 60 / IOPROG_60.ISO / soft / c++ / gsl-1.1.1-setup.exe / {app} / src / integration / qc25f.c < prev    next >
Encoding:
C/C++ Source or Header  |  2001-04-18  |  4.2 KB  |  159 lines
  1. /* integration/qc25f.c
  2.  * 
  3.  * Copyright (C) 1996, 1997, 1998, 1999, 2000 Brian Gough
  4.  * 
  5.  * This program is free software; you can redistribute it and/or modify
  6.  * it under the terms of the GNU General Public License as published by
  7.  * the Free Software Foundation; either version 2 of the License, or (at
  8.  * your option) any later version.
  9.  * 
  10.  * This program is distributed in the hope that it will be useful, but
  11.  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12.  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13.  * General Public License for more details.
  14.  * 
  15.  * You should have received a copy of the GNU General Public License
  16.  * along with this program; if not, write to the Free Software
  17.  * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  18.  */
  19.  
  20. struct fn_fourier_params
  21. {
  22.   gsl_function *function;
  23.   double omega;
  24. };
  25.  
  26. static double fn_sin (double t, void *params);
  27. static double fn_cos (double t, void *params);
  28.  
  29. static void
  30. qc25f (gsl_function * f, double a, double b, 
  31.        gsl_integration_qawo_table * wf, size_t level,
  32.        double *result, double *abserr, double *resabs, double *resasc);
  33.  
  34. static void
  35. qc25f (gsl_function * f, double a, double b, 
  36.        gsl_integration_qawo_table * wf, size_t level,
  37.        double *result, double *abserr, double *resabs, double *resasc)
  38. {
  39.   const double center = 0.5 * (a + b);
  40.   const double half_length = 0.5 * (b - a);
  41.   const double omega = wf->omega ;
  42.   
  43.   const double par = omega * half_length;
  44.  
  45.   if (fabs (par) < 2)
  46.     {
  47.       gsl_function weighted_function;
  48.       struct fn_fourier_params fn_params;
  49.  
  50.       fn_params.function = f;
  51.       fn_params.omega = omega;
  52.  
  53.       if (wf->sine == GSL_INTEG_SINE) 
  54.     {
  55.       weighted_function.function = &fn_sin;
  56.     }
  57.       else
  58.     {
  59.       weighted_function.function = &fn_cos;
  60.     }
  61.  
  62.       weighted_function.params = &fn_params;
  63.  
  64.       gsl_integration_qk15 (&weighted_function, a, b, result, abserr,
  65.                 resabs, resasc);
  66.       
  67.       return;
  68.     }
  69.   else
  70.     {
  71.       double *moment;
  72.       double cheb12[13], cheb24[25];
  73.       double result_abs, res12_cos, res12_sin, res24_cos, res24_sin;
  74.       double est_cos, est_sin;
  75.       double c, s;
  76.       size_t i;
  77.  
  78.       gsl_integration_qcheb (f, a, b, cheb12, cheb24);
  79.  
  80.       if (level >= wf->n)
  81.     {
  82.           /* table overflow should not happen, check before calling */
  83.           GSL_ERROR_VOID("table overflow in internal function", GSL_ESANITY);
  84.     }
  85.  
  86.       /* obtain moments from the table */
  87.  
  88.       moment = wf->chebmo + 25 * level;
  89.  
  90.       res12_cos = cheb12[12] * moment[12];
  91.       res12_sin = 0 ;
  92.  
  93.       for (i = 0; i < 6; i++)
  94.     {
  95.       size_t k = 10 - 2 * i;
  96.       res12_cos += cheb12[k] * moment[k];
  97.       res12_sin += cheb12[k + 1] * moment[k + 1];
  98.     }
  99.  
  100.       res24_cos = cheb24[24] * moment[24];
  101.       res24_sin = 0 ;
  102.  
  103.       result_abs = fabs(cheb24[24]) ;
  104.  
  105.       for (i = 0; i < 12; i++)
  106.     {
  107.       size_t k = 22 - 2 * i;
  108.       res24_cos += cheb24[k] * moment[k];
  109.       res24_sin += cheb24[k + 1] * moment[k + 1];
  110.       result_abs += fabs(cheb24[k]) + fabs(cheb24[k+1]);
  111.     }
  112.  
  113.       est_cos = fabs(res24_cos - res12_cos);
  114.       est_sin = fabs(res24_sin - res12_sin);
  115.  
  116.       c = half_length * cos(center * omega);
  117.       s = half_length * sin(center * omega);
  118.  
  119.       if (wf->sine == GSL_INTEG_SINE)
  120.     {
  121.       *result = c * res24_sin + s * res24_cos;
  122.       *abserr = fabs(c * est_sin) + fabs(s * est_cos);
  123.     }
  124.       else
  125.     {
  126.             *result = c * res24_cos - s * res24_sin;
  127.       *abserr = fabs(c * est_cos) + fabs(s * est_sin);
  128.     }
  129.       
  130.       *resabs = result_abs * half_length;
  131.       *resasc = GSL_DBL_MAX;
  132.  
  133.       return;
  134.     }
  135. }
  136.  
  137. static double
  138. fn_sin (double x, void *params)
  139. {
  140.   struct fn_fourier_params *p = (struct fn_fourier_params *) params;
  141.   gsl_function *f = p->function;
  142.   double w = p->omega;
  143.   double wx = w * x;
  144.   double sinwx = sin(wx) ;
  145.   return GSL_FN_EVAL (f, x) * sinwx;
  146. }
  147.  
  148. static double
  149. fn_cos (double x, void *params)
  150. {
  151.   struct fn_fourier_params *p = (struct fn_fourier_params *) params;
  152.   gsl_function *f = p->function;
  153.   double w = p->omega;
  154.   double wx = w * x;
  155.   double coswx = cos(wx) ;
  156.   return GSL_FN_EVAL (f, x) * coswx ;
  157. }
  158.  
  159.